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Extensive chromosomal reshuffling drives evolution 
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Sexual recombination drives genetic diversity in euicaryotic genomes and fosters adaptation to novel environmental 
challenges. Although strictly asexual microorganisms are often considered as evolutionary dead ends, they comprise 
many devastating plant pathogens. Presently, it remains unknown how such asexual pathogens generate the genetic 
variation that is required for quick adaptation and evolution in the arms race with their hosts. Here, we show that 
extensive chromosomal rearrangements in the strictly asexual plant pathogenic fungus Verticillium dahliae establish highly 
dynamic lineage-specific [LS) genomic regions that act as a source for genetic variation to mediate aggressiveness. We show 
that such LS regions are greatly enriched for in planta-expressed effector genes encoding secreted proteins that enable host 
colonization. The LS regions occur at the flanks of chromosomal breakpoints and are enriched for retrotransposons and 
other repetitive sequence elements. Our results suggest that asexual pathogens may evolve by prompting chromosomal 
rearrangements, enabling rapid development of novel effector genes. Likely, chromosomal reshuffling can act as a general 
mechanism for adaptation in asexually propagating organisms. 



[Supplemental material is available for this article.] 

Sexual reproduction is thought to enable responses to environ- 
mental challenges in nature and fosters evolution of pathogen 
genomes (Williams 1975; Maynard Smith 1978; Colegrave 2002; 
Goddard et al. 2005; Heitman 2006; de Visser and Elena 2007). 
Consequently, pathogen populations that undergo regular sexual 
reproduction are thought to pose a great risk to agriculture because 
they can recombine alleles that contribute to virulence in the face 
of dynamic environmental conditions (McDonald and Linde 
2002). However, sexual reproduction comes at a cost because two 
compatible individuals need to locate each other to generate off- 
spring, and fitness decreases due to break up of coadapted com- 
binations of interacting alleles may occur (Agrawal 2006; Heitman 
2006; Heitman et al. 2007; de Visser and Elena 2007). Sexual 
compatibility between individuals depends on so-called mating 
types in fungi. In heterothallic fungi, mating is restricted to ge- 
netically different individuals of opposite mating t3^es, whereas in 
homothallic fungi, mating and sexual reproduction can occur 
between all individuals, including itself (Billiard et al. 2012). Al- 
though most fungi are able to undergo asexual and sexual re- 
production, —20% of all fungal phyla reproduce strictly asexually 
(Heitman et al. 2007). Such strictly asexual organisms are thought 
to be less flexible than sexual ones, relying solely on random 
mutations to adapt to changing environments. Asexual organisms 
are often considered as evolutionary dead ends (Burt 2000; McDonald 
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and Linde 2002). This is mainly due to the absence of meiotic re- 
combination resulting in increased accumulation of deleterious 
mutations, an effect known as Muller's ratchet (Felsenstein 1974), 
and a decreased ability to adapt to changing environmental con- 
ditions by generation of novel genetic combinations (Heitman 
2006). Especially for pathogenic species that are continuously in- 
volved in an arms race with their hosts that try to ward off invading 
pathogens, quick adaptation in order to coevolve with the host 
immune system is essential for evolutionary success (Raffaele and 
Kamoun 2012). Nonetheless, many destructive plant pathogenic 
fungi are known as strictly asexually reproducing organisms, such as 
the vascular wilt fungus Verticillium dahliae (Fradin and Thomma 
2006). 

V. dahliae is a soil-borne broad host-range plant pathogen that 
invades the water-conducting xylem vessels of susceptible plant 
species to cause vascular wilt disease (Fradin and Thomma 2006; 
Klosterman et al. 2009). Hundreds of dicotyledonous plant species 
can be infected by V. dahliae, including many economically im- 
portant crops such as lettuce, cotton, and tomato (Bolek et al. 
2005; Fradin and Thomma 2006; Klosterman et al. 2009; Atallah 
et al. 2010). Although two mating types have been found in 
V. dahliae, it is predicted that propagation is solely asexual (clonal) 
in most if not all populations, as a sexual cycle has never been 
observed and the ratio between both mating types is greatly 
skewed toward one of the two (Usami et al. 2008, 2009; Atallah 
et al. 2010; Inderbitzin et al. 2011). Only in a few plant species has 
monogenic resistance toward Verticillium wilt been described. 
These include tomato, for which a Verticillium resistance locus 
(named Ve) has been identified that mediates resistance against 
race 1 strains of V. dahliae and V. albo-atrum, and that thus become 
avirulent (unable to cause disease). Strains that are not contained 
by this Ve locus remain virulent (able to cause disease) and are 
assigned to race 2 (Schaible et al. 1951; Kawchuk et al. 2001; Fradin 
et al. 2009). The resistance that is established by the Ve locus is 
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mediated by the Vel gene that encodes a receptor-like protein-type 
cell surface receptor (Fradin et al. 2009) that activates plant immu- 
nity upon perception of a corresponding pathogen-derived ligand. 
Recently through comparative genomics the race 1 -specific effector 
Avel (for Avirulence on Vel tomato) was identified that activates 
Vel-mediated resistance in tomato (de Jonge et al. 2012). Remarkably 
the Avel effector contributes to pathogen aggressiveness on tomato 
plants that lack Vel, demonstrating that Avel acts as a virulence 
factor of V. dahliae (de Jonge et al. 2012). Race 2 strains lack the Avel 
gene, and are consequently less aggressive on tomato plants that 
lack Vel when compared with race 1 strains (Amen and Shoemaker 
1985; Paternotte and van Kesteren 1993; de Jonge et al. 2012). 

In this study, we compared the genomes of 11 recently se- 
quenced V. dahliae strains (Klosterman et al. 2011; de Jonge et al. 
2012) that have been isolated from various geographical locations 
and hosts, and that are all, except for one, pathogenic on tomato. 
By assessing genetic diversity within this population we aimed to 
study the impact of asexual reproduction on the evolution of host 
adaptation and virulence in this plant pathogenic fungus. 

Results 

Identification of the core genome reveals extensive 
chromosomal rearrangements 

With comparative genomics, we examined genetic diversity in a 
population of nine recently sequenced tomato pathogenic strains 
and one nonpathogenic strain of the asexual fungus Verticillium 
dahliae (Supplemental Table 1; de Jonge et al. 2012) together with 
the genome of V dahliae strain VdLsl7 as a reference (Klosterman 
et al. 201 1). Read coverage mapping in 1-kb windows of all sequenced 
V. dahliae strains over the VdLsl7 reference genome revealed a 
"core genome" that is shared by all strains, encompassing —29 Mb 
of sequence containing 9471 genes. The amount of single nucle- 
otide polymorphisms (SNPs) within the core genome ranged from 
5445 (JR2; -99.98% identity) to 163,602 (St. 100; -99.5% identity) 
SNPs per strain (Table 1), collectively amounting to 236,785 non- 
redundant polymorphic sites. Of these, 78,342 (32.9%) occurred in 
protein-coding regions, of which 55% were synonymous, not af- 
fecting the protein sequence, and 45% were nonsynonymous 
(Table 1). To determine selection strength, the ratio of non- 
synonymous (d^) substitutions per nonsynonymous site (Xa) to 
the number of synonymous substitutions (ds) per synonymous site 
(Xs) was calculated within the coding regions for each of the 9471 
core genes (Stukenbrock et al. 2010). A surprisingly low number of 
29 genes was found to be under positive selection (KJK^ > 1; P < 
0.01). About half of these (13) encode proteins with a known 



function, including cellulose binding and regulation of transcrip- 
tion, whereas the other half encodes proteins of unknown func- 
tion, of which only three encode secreted proteins that are can- 
didate effectors (Supplemental Table 2; de Jonge et al. 2011). 

To infer evolutionary relations within the population, all 
236,785 SNP positions were used to construct a phylogenetic tree 
(Fig. 1). Despite the high degree of conservation between the 
strains that were included in the population (maximum nucleotide 
divergence within the population is only 0.5%), strain JR2 formed 
a clearly separate cluster with the VdLsl 7 reference strain, differing 
only by one SNP in -6000 bp (0.02%). Despite the high degree of 
identity (99.98%), the two strains show differential aggressiveness 
on tomato. Strain JR2 belongs to race 1 and shows a high degree of 
aggressiveness on tomato lacking Vel, while VdLsl 7 belongs to 
race 2 and displays a relatively mild aggressiveness on tomato 
lacking Vel . These two strains were analyzed in further detail. To 
improve the de novo JR2 genome assembly, mate-pair sequencing 
of a 5-kb insert library was performed, leading to a drastically im- 
proved contig N50 from 59.4 kb to 1.0 Mb, and a decreased 
number of contigs from 4753 to 267. We subsequently used optical 
mapping to obtain an accurate whole-genome assembly and 
placed —30% of the contigs on eight scaffolds, covering —94% of 
the total assembly (34 Mb) (Supplemental Table 3). While the 
chromosome lengths of the reference strain were found to vary 
between 3.1 and 6.0 Mb (Klosterman et al. 2011), JR2 chromosome 
sizes differ and range up to 8.9 Mb (Supplemental Table 4). Thus, 
despite the observation that the two strains are 99.98% identical in 
all genomic regions that could be aligned, the karyotype of these two 
strains is completely different. This observation urged us to in- 
vestigate collinearity between the chromosomes of the two strains. 

Unanticipated, pairwise alignments to identify collinear 
(synteny) blocks between the two strains revealed repeated in- 
terruptions of syntenic regions by intra- and interchromosomal 
rearrangements (Fig. 2). In total, 11 intra- and 17 interchro- 
mosomal rearrangements were identified between the two strains 
(Fig. 2A,B). While most of these breakpoints were associated with 
assembly gaps due to the occurrence of repeats in the flanking 
genomic areas (Fig. 2B,C; Supplemental Table 6), three breakpoints 
could be pinpointed to the nucleotide exactly, and thus could be 
experimentally confirmed with PGR (Fig. 1; Supplemental Fig. 1; 
Supplemental Tables 5, 6). We subsequently extended the analysis 
of chromosomal rearrangements to all sequenced strains by 
screening for breakpoints within alignments of the assembled 
contigs for each of the strains to the reference strain. Although 
identification of such syntenic breakpoints in the small insert li- 
brary genome assemblies is hampered by the short contig lengths, 
three to eight synteny breakpoints could be identified in silico for 



Table 1. Summary of SNPs when compared with V. dahliae reference strain VdLsl 7 

Strain #SNPs #Unique SNPs #SNPs in intergenic region #SNPs in introns #SNPs in exons c/n^ ds^ 



JR2 


5,445 


947 


3,193 


636 


1,563 


1,092 


444 


CBS381 .66 


117,364 


166 


67,640 


8,775 


39,062 


16,702 


21,974 


St14.01 


117,704 


274 


67,894 


8,786 


39,161 


16,744 


22,030 


St. 100 


163,602 


81061 


94,516 


12,108 


54,219 


23,350 


30,361 


DVD-3 


118,528 


211 


68,773 


8,716 


39,135 


16,662 


22,098 


DVD-31 


119,025 


316 


68,973 


8,786 


39,381 


16,827 


22,170 


DVD-1 61 


117,277 


178 


67,831 


8,656 


38,776 


16,526 


21,879 


DVD-S26 


117,307 


235 


67,653 


8,764 


39,079 


16,743 


21,947 


DVD-S29 


122,057 


5552 


70,801 


9,204 


40,219 


17,175 


22,647 


DVD-S94 


1 1 9,060 


205 


68,993 


8,782 


39,303 


16,763 


22,161 


All 


236,785 


89,145 


136,705 


1 7,421 


78,342 


34,007 


43,563 



^Number of nonsynonymous (cfn) and synonymous substitutions (ds). 
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Figure 1. Population structure and chromosome karyotype of sequenced Verticillium dahliae isolates. An unrooted maximum likelihood phylogenetic 
tree based on concatenation of 236,785 SNP sites relative to reference strain VdLsl 7 is shown. Evolutionary distances based on the Jukes-Cantor method 
and bootstrap support (%) are indicated at the nodes. Race 2 strains are indicated by asterisks. Pulsed-field gel electrophoresis of sequenced V. dahliae 
strains are shown as well as the results of PCR analyses on structural rearrangements SR2, SR3, and SR4 using JR2-specific primer sets (J2, J3, J4) and VdLsl 7- 
specific primer sets (L3, L4), highlighting significant chromosome length and structure polymorphisms between isolates. Chromosomal DNA of Schizo- 
saccharomyces pombe and Hansenula wingei were loaded as size makers. 



each of these strains when compared with the VdLsl 7 reference 
genome, which roughly occur in the regions of the previously 
identified breakpoints between VdLsl 7 and JR2. The occurrence 
of several of these breakpoints was confirmed with PCR (Fig. 1; 
Supplemental Fig. 1). Moreover, pulsed-field gel electrophoresis 
confirmed considerable chromosome length polymorphism be- 
tween all strains (Fig. 1) and demonstrates that intra- and in- 
terchromosomal rearrangements frequently occur within the ge- 
nomes of V. dahliae strains despite their high degree of sequence 
conservation. 

Particularly retrotransposons have been implicated in genome 
rearrangements through homologous recombination between ele- 
ments or by causing chromosomal breaks during excision or in- 
sertion (Mieczkowski et al. 2006; Maxwell et al. 2011). To investigate 
a potential involvement of retrotransposons in the genomic rear- 
rangements of V. dahliae, we identified all repetitive elements in the 
genome sequences of strains VdLsl 7 and JR2, and found that both 
genomes contain only a low-repeat content of —4% (Klosterman 
et al. 201 1; Supplemental Table 7). Of the different classes of repeats, 
we only observed significant correlation (P < 0.01) between synteny 
breakpoints and long terminal repeat (LTR) retrotransposons (Sup- 
plemental Tables 8-11). 

The pan-genome of V. dahliae: Lineage-specific regions 

In addition to the core —30 Mb genome, all strains carried up to ~4 
Mb of genome sequence that was unique or shared by only a subset 
of strains, composing a highly dynamic lineage-specific (LS) region 



of the genome encoding up to 1000 genes (Table 2; Klosterman 
et al. 2011), which is enriched for various types of transposable 
elements (Fig. 2B,C; Klosterman et al. 2011; Amyotte et al. 2012). 
These LS genomic regions are correlated (P < 0.01) with syntenic 
breakpoints in VdLsl 7 and JR2 (Supplemental Tables 6, 12). In 
VdLsl 7 as well as JR2 these regions are typically devoid of house- 
keeping genes that are required for the maintenance of basic cel- 
lular functions such as transcription and translation (Klosterman 
et al. 2011). Intriguingly, th^Avel gene is located within the ~4-Mb 
LS region of race 1 strains, suggesting that these regions contribute 
to niche adaptation, i.e., pathogenicity on plant hosts (de Jonge 
et al. 2012). 

To further establish the role of LS regions in V. dahliae path- 
ogenicity, we mined these regions, as well as the core genomic 
regions in VdLsl 7 and JR2, for candidate effectors (i.e., small, se- 
creted proteins of unknown function). Typically, effector genes 
encode small proteins that are characterized by the presence of an 
N-terminal signal peptide that can be predicted computationally 
(de Jonge 2012). However, we did not observe a remarkable dif- 
ference in the number of small secreted proteins between the core 
and LS genomic regions (Table 2). Similarly, no remarkable differ- 
ence was found in the number of protein domains classified in the 
Pfam database (Finn et al. 2008). In a number of plant pathogenic 
oomycetes it was observed that effector genes are found in repeat- 
rich, gene-sparse regions that can be visualized by plotting the 
length of the 5'- and 3 '-flanking intergenic regions (Raffaele et al. 
2010). In the case of V. dahliae, we did not find a significant difference 
in gene density between core and LS genomic regions (Table 2; Sup- 
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Alignment to unpositioned scaffolds (Un) 



Figure 2. Whole-genome alignment of Verticillium dahliae strains VdLsl 7 and JR2 reveals extensive chromosomal rearrangements. (A) Whole-genome 
dot-plot comparison with forward-forward alignments (black) and inversions (blue). Red triangles mark syntenic breakpoints. (Un) Unplaced contigs 
during optical mapping. (B) Global view of synteny alignments between V. dahliae strains VdLsl 7 and JR2 and the distribution of LTR retrotransponsons 
and lineage-specific regions. VdLsl 7 chromosomes are shown as reference. For each chromosome, row I represents genomic scaffolds (black) on the 
chromosome separated by scaffold breaks (gray); row II displays syntenic alignment of JR2 chromosomes, highlighting 28 major intra- and in- 
terchromosomal rearrangements that are marked by red (inter) and blue (intra) arrows above each breakpoint; row III represents the density of LTR 
retrotransposable elements calculated in 1 0-kb windows; and row IV depicts lineage-specific regions that are found in only a subset of strains. (C) Circos 
diagram illustrating collinear blocks with alignments between VdLsl 7 (gray) and JR2 (blue) chromosomes, sequence gaps, sequences aligning to 
unpositioned scaffolds (yellow), lineage-specific sequences (red), repeat density (% coverage of 1 0-kb window), and GC % (per 1 0-kb window). 



plemental Fig. 2). However, when assessing transcription (RNA-seq) 
data of V. dahliae-infected N. benthamiana (de Jonge et al. 2012; 
Faino et al. 2012), we identified 19 V. dahliae genes (P < 0.001) that 
were induced to the same level, or higher, as the previously char- 
acterized ylvei effector gene [logio(FC) > 2.15; 8-d post-inoculation] 



when compared with expression in in vitro-cultured mycelium. 
Remarkably, eight of the 19 genes are located within the ~3-Mb LS 
region, while the remaining 11 genes reside in the ~30-Mb core 
genome, identifying a significant overrepresentation of in planta- 
induced genes in the LS regions (Fig. 3). Of the eight highly 



Table 2. Summary of genome size and gene numbers for each sequenced strain, split by the core and lineage-specific (LS) genome 

Genome size (Mb) Number of genes (%) Secretome 



Strain Total Excl. gaps Repeats Core LS Total Secreted Core LS Core/secreted LS/secreted Core LS 



St.14.01 


34.3 


33.9 


2.3 


30.4 


3.4 


9,854 


979 


9,341 


513 


941 


38 


10.1 


7.4 


DVD-3 


33.8 


33.3 


1.9 


30.5 


2.8 


9,881 


972 


9,335 


546 


937 


35 


10.0 


6.4 


CBS381 .66 


34.5 


34.2 


2.6 


30.4 


3.8 


9,834 


1,005 


9,230 


604 


961 


44 


10.4 


7.3 


St. 100 


35.0 


34.6 


0.8 


30.1 


4.5 


9,822 


976 


9,193 


629 


916 


60 


10.0 


9.5 


JR2 


37.5 


33.5 


1.5 


30.1 


3.4 


10,985 


1,085 


10,098 


887 


1,005 


80 


10.0 


9.0 


DVD-S94 


34.6 


34.2 


2.6 


30.5 


3.8 


9,846 


989 


9,289 


557 


951 


38 


10.2 


6.8 


DVD-31 


33.8 


33.4 


2.4 


30.4 


3.0 


9,867 


985 


9,505 


362 


951 


34 


10.0 


9.4 


DVD-S26 


35.0 


34.6 


2.3 


30.4 


4.2 


10,011 


970 


9,336 


675 


930 


40 


10.0 


5.9 


DVD-S29 


33.7 


33.1 


2.0 


30.2 


2.9 


9,773 


969 


9,275 


498 


939 


30 


10.1 


6.0 


DVD-1 61 


33.8 


33.4 


2.2 


30.4 


3.0 


9,852 


997 


9,427 


425 


967 


30 


10.3 


7.1 


VdLsl 7 


36.9 


33.0 


1.3 


29.2 


3.7 


10,535 


1,055 


9,530 


1,005 


974 


81 


10.2 


8.1 



1274 Genome Research 



www.genome.org 



Verticillium chromosomal reshuffling 




Variable 
(887 genes) 





In planta induced genes 
logio(FC)>2 (19 genes) 



XLOC_009059 

XLOC_008951 

evm. model. contiglBGB. 344 

evm. model. contigUn. 113 

evm. model. contig34339. 933 

XLOC_001220 

evm. model. contig34339. 698 

XLOC_00170 

evm. model. contigUn. 128 

evm. model. contig4566. 2197 

XLOC_001183 

evm. model. contig44686.971 
evm. model. contig44686.970 
evm.model.contigl569.601 
evm. model. contig34339. 864 
evm. model. contig44576.80 
evm. model, con tig45503.424 
evm. model. contig4566. 2834 
Avel 



CBS381.6( 
DVD-S26 


DVD-S29 


DVD-161 


VdLsl7 


DVD-3 


St.14.01 


St.lOO 


DVD-S94 


DVD-31 


9 












53/72 












■ 












































37/56 
































1 











































































































































12 16 



DPI 

Figure 3. Lineage-specific genomic regions of Verticillium dahliae strain JR2 are enriched for in planta-expressed genes. Deep transcriptome sequencing 
of infected Nicotiana benthamiana plants harvested between 4- and 1 6-d post-inoculation (DPI), and fungus cultured on Czapek Dox medium (C) led to 
the identification of 19 highly in planta-induced genes [logio(FC) > 2.1 5]. Of these, eight are located in lineage-specific (LS) regions (in blue), while 1 1 are 
in the core genome (in red), demonstrating significant overrepresentation of in planta-expressed genes in LS regions. The table shows the presence (black) 
or absence (white) of in planta-induced JR2 genes within the set of sequenced strains and indicates hypothetical and putatively secreted proteins. Numbers 
in the two gray cells indicate amino acid identity and conservation, respectively, for two divergent homologs. 



induced genes that reside in LS regions, five (including Avel) en- 
code putative secreted effectors (Fig. 3). Moreover, considering a 
logio(FC) > 3 or logio(FC) > 4, these ratios are even more skev^^ed v^^ith 
six out of 14 and five out of six in planta-induced genes residing in 
LS regions, respectively. 

Lineage-specific genomic regions determine V. dahliae virulence 
and niche adaptation 

To assess the biological significance of the significant over- 
representation of in planta-induced genes in the LS regions, we 
pursued targeted deletion of tv^^o randomly chosen effector genes 
of the four effector gene candidates in addition to Avel that are most 
highly expressed in planta; XLOC_009059 and XLOC_008951. 
XLOC_009059 does not occur in any of the other strains except for 
St.lOO, v^hich contains a highly diversified copy (53% and 72% 
amino acid identity and conservation, respectively, as determined 
by BLASTp) (Fig. 3) that is unlikely to encode a functional protein 
due to the high number of predicted stop codons. On the other 
hand, XLOC_008951 was found in four of the 11 strains that v^^ere 
sequenced (Fig. 3). Remarkably, deletion of each of them resulted 
in significantly reduced disease development, confirming that 
they encode bona fide effectors required for virulence on tomato 
(Fig. 4A,B). Thus, we conclude that the LS regions of strain JR2 
significantly contribute to virulence. 



Similar to JR2, strain VdLsl7 also contains a ~4-Mb LS region 
(Table 2; Klosterman et al. 2011) that was assessed for unique ef- 
fector genes crucial for pathogenicity. All V. dahliae strains that 
weie sequenced contain six genes that encode LysM effectors v\^ithin 
their core genome {VDAG_00902] VDAG_03096; VDAG_04781; 
VDAG_06426; VDAG_06998; VDAG_08171). LysM effectors have 
previously been implicated in sequestration of chitin fragments 
to prevent these fragments from stimulating host immunity (de 
Jonge and Thomma 2009; de Jonge et al. 2010; Marshall et al. 201 1; 
Mentlak et al. 2012). VdLslZ contains an additional LysM effector 
gene {VDAG_05180) within a LS region, encoding a LysM effector 
that contains two LysM domains. Remarkably, of the seven LysM 
effectors in the VdLsl7 genome, only this LysM effector gene is 
significantly expressed in planta (Fig. 4C). Furthermore, targeted 
deletion of this effector revealed that it is indeed required for vir- 
ulence of the VdLslZ strain (Fig. 4D,E). Thus, like in JR2, the LS 
regions of strain VdLsl7 contribute to virulence. Notably, like the 
race-l-specific ettectoiAvel (de Jonge et al. 2012), VDAG_05180 as 
well as V. dahliae ]R2 effector genes XLOC_009059, XLOC_008951 
and the candidate effector genes XLOC_00170 and evm. model. 
contigUn. 128 are flanked by repetitive elements such as LTR retro- 
transposons and FotI transposases (Fig. 5; Supplemental Figure 3), 
suggesting that the presence/absence polymorphism of all 
of these effector genes in V. dahliae is mediated by the flexibility 
and instability of the genomic context as a consequence of 
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WT (VdLs17) 



AVDAG_05180 



Figure 4. Lineage-specific genomic regions of Verticillium dahliae harbor effectors that are required 
for virulence on tomato. (A) Two independent XLOC_009059 deletion strains (AXLOC_009059) show 
compromised virulence on tomato (cv. MoneyMaker), evidenced by reduced stunting when compared 
with inoculation with wild-type V. dahliae (WT, JR2), increased canopy area, and reduced V. dahliae 
biomass accumulation. (B) Two independent XLOC_008951 deletion strains (AXLOC_008951) show 
compromised virulence on tomato (cv. MoneyMaker), evidenced by reduced stunting when compared 
with inoculation with wild-type V. dahliae (WT, JR2), increased canopy area, and reduced V. dahliae 
biomass accumulation. (C) Two independent VDAG_051 80 deletion strains (AVDAG_051 80) show 
compromised virulence on tomato (cv. Motelle), evidenced by reduced stunting when compared with 
inoculation with wild-type V. dahliae (WT, VdLsl 7). Photographs are taken at 1 2 DPI. (D) Expression of 
VDAG_051 80\r\ V. cfoMoe strain VdLsl / during infection of N/cof/onobenf/?o/?7/ono between 4- and 16-d 
post-inoculation (DPI). (£) Reduced V. dahliae biomass in plants inoculated with two independent 
VDAG_05180 deletion strains when compared with wild-type V. dahliae (WT, VdLsl 7) at 8 DPI. Error 
bars represent standard error of three replicate experiments. 



these repetitive elements (Klosterman et al. 2011; Amyotte et al. 
2012). 

Discussion 

So far, it remains unknown how asexually propagating organisms, 
including many fungal species, are able to adapt to changing en- 
vironmental conditions. The lack of sexual reproduction, and 
consequently meiotic recombination, is thought to greatly affect 
the potential to generate novel genetic combinations. Especially, 
pathogenic microbes involved in an intimate interaction with 
a host need to continuously adapt their effector gene repertoire to 



counteract host defense responses that 
have evolved in the evolutionary arms 
race with the pathogen. Here, we exam- 
ined the genetic diversity in a population 
of the plant pathogenic fungus V. dahliae 
that propagates solely through asexual 
reproduction. As expected for an asexually 
reproducing organism, we only found a 
low degree of genome-wide nucleotide 
diversity between V. dahliae strains in the 
population, ranging between 0.02% and 
0.5%. For comparison, in a study on the 
sexually reproducing fungal wheat path- 
ogen Zymoseptoria tritici (synonym of 
Mycosphaerella graminicola), the average 
genome-wide nucleotide diversity was 
determined to range between 2% and 4% 
(Stukenbrock et al. 201 1). Despite the low 
degree of genetic diversity, we identified 
29 genes that are under positive selection 
within the V. dahliae population (Sup- 
plemental Table 2), suggesting that these 
genes are under adaptive evolution. This 
number is comparable to the number of 
genes under positive selection in Z. tritici 
when compared with its wild sister spe- 
cies (Stukenbrock et al. 2011). However, 
only three of these V. dahliae genes en- 
code secreted proteins that are candidate 
effectors, suggesting that molecular evo- 
lution at the nucleotide level does not 
play a major role in host adaptation of 
this asexual pathogen. 

Various mechanisms have been de- 
scribed that facilitate rapid development 
of novel effector genes in pathogenic 
microbes, including diversity at genomic 
locations enriched for transposons, mu- 
tation, and recombination in subtelomeric 
regions (McDonagh et al. 2008; Chuma 
et al. 2011), coregulated gene clusters 
(Pallmer and Keller 2010; Schirawski et al. 
2010), small dispensable chromosomes 
(Coleman et al. 2009; Ma et al. 2010; 
Stukenbrock et al. 2010; Goodwin et al. 
2011; Raffaele and Kamoun 2012), gene 
sparse regions (Raffaele et al. 2010), AT- 
rich isochore-like regions (van der Wouw 
et al. 2010; Rouxel et al. 2011), genome 
hybridization (Stukenbrock et al. 2012), 
and horizontal gene transfer (Friesen et al. 2006; de Jonge et al. 
2012; Gardiner et al. 2012). However, most of these mechanisms 
have been described in species that can reproduce sexually, and of 
which the genomes were shaped by repeat-driven expansion 
(Raffaele and Kamoun 2012). 

Here we describe the observation of a novel mode of evolu- 
tion of pathogenicity in the asexual plant pathogenic fungus 
V. dahliae, of which only 4% of the genome is repetitive (Supple- 
mental Table 7; Klosterman et al. 2011; Amyotte et al. 2012). This is 
significantly less when compared with the genomes of other fila- 
mentous pathogens, such as the fungi Magnaporthe oryzae (10%), 
Z. tritici (17%), Fusarium oxysporum (28%), Cladosporium fulvum 
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Figure 5. Genomic context of the lineage-specific Verticillium dahliae LysM effector VDAG_051 80. (A) Genomic location of VDAG_051 80 in V. dahliae 
strain VdLsl 7 revealing the presence of flanking gaps, repeats, and a predicted overlapping, but inactive LTR retrotransposon. RepeatModelerfamilies rnd- 
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(47%), Blumeria graminis (64%), and the oomycetes Hyalo- 
peronospora arabidopsidis (42%) and Phytophthora infestans (74%) 
(Dean et al. 2005; Haas et al. 2009; Baxter et al. 2010; Ma et al. 
2010; Spanu et al. 2010; Goodwin et al. 2011; de Wit et al. 2012). 
Comparative analyses between highly similar V. dahliae strains 
revealed numerous intra- and interchromosomal rearrangements 
and extensive karyotype variation. This is remarkable, as it is 
generally assumed that genome rearrangement increases with in- 
creasing sequence divergence. Previous analyses of fungal ge- 
nomes have revealed that although large syntenic regions with 
orthologous genes arranged in the same orientation across species 
(known as macrosynteny) are rare, chromosomal gene content 
is typically preserved among Ascomycetes, particularly in the 
Dothidiomycete class (Goodwin et al. 201 1; Hane et al. 201 1). Such 
preservation of chromosomal gene content without conservation 
of gene order or orientation is referred to as mesosynteny, and is 
characterized by the occurrence of numerous intrachromosomal 
rearrangements, presumably arising as the consequence of fre- 
quent inversions during meiosis rather than interchromosomal 
rearrangements (Hane et al. 2011). The extensive number of in- 
terchromosomal rearrangements between V. dahliae isolates in- 
dicates that mesosynteny is not crucial in V. dahliae. 

Unlike previous observations in a number of fungal plant 
pathogens, including F. oxysporum, N. haematococca, and Z. tritici, 
which display karyotype variation as a consequence of the differ- 
ential presence of small dispensable chromosomes (Coleman et al. 
2009; Ma et al. 2010; Stukenbrock et al. 2010; Goodwin et al. 201 1; 
Raffaele and Kamoun 2012), the karyotype variation that we ob- 
served between individual isolates of V. dahliae concerns chromo- 
some length polymorphisms, and no indication for the occurrence 
of small dispensable chromosomes is found in V. dahliae. The ob- 
served variation in this study and previous ones among fungal 
karyotypes (Zolan 1995) is the result of complex chromosomal 
rearrangements. It has previously been suggested that karyotype 
variation, rather than being a mechanism of adaptation to gener- 
ate novel virulence traits, occurs because nondeleterious genomic 
rearrangements are maintained due to the absence or rarity of 
a sexual cycle (Talbot et al. 1993; Zolan 1995). Our data challenge 
this hypothesis by showing that chromosomal plasticity, evi- 
denced by extensive targeted chromosomal rearrangements and 
karyotype variability, is genetically non-neutral as it induces local 
sequence variation at syntenic breakpoints and increases adaptive 
capability, as demonstrated by the significant enrichment for 
novel effector genes and acquisition of the race 1 -specific effector 
Avel, possibly through horizontal gene transfer and subsequent 
loss in V. dahliae race 2 isolates (de Jonge et al. 2012). Moreover, our 
data reinforces the notion that V. dahliae strictly relies on asexual 
reproduction, as extensive karyotype variation prevents correct 
pairing of homologous chromosomes during meiosis (Kistler and 
Miao 1992). 

Our phylogenetic analysis based on whole-genome sequence 
data (Fig. 1) has demonstrated that race 1 and race 2 isolates of 
V. dahliae do not form separate clades, and the fact that we did not 
observe any sequence variation for Avel within and between Ver- 
ticillium species (de Jonge et al. 2012) strongly suggests that Avel 
has been acquired (in race 1) or lost (in race 2) several times. Sim- 
ilarly, we identified additional effector genes within the plastic 
genomic regions that are unique or shared by only a subset of 
strains and like Avel; these are conserved and not monophyletic, 
suggesting that they are frequently lost and/or readily exchanged 
between individual strains as has been proposed for effector genes 
of the Avr-Pita family in M. oryzae (Chuma et al. 201 1). We speculate 



that frequent presence/absence polymorphisms of effector genes 
in these plastic regions is mediated by the flexibility and instability 
of these regions, reflected by the enrichment for transposable ele- 
ments and association with chromosomal rearrangements (Amyotte 
et al. 2012). Moreover, we note that particular transposable ele- 
ments flank a number of the in planta highly induced effector 
genes located in the LS regions and speculate that, besides pro- 
viding genetic flexibility, these elements might also affect the level 
of gene expression. An association between in planta-expressed 
effector genes and a particular transposable element was recently 
found for the tomato pathogenic vascular wilt fungus Fusarium 
oxysporum f. sp. lycopersici, in which a repetitive miniature Impala 
(mimp) element was found to be present in the promoter region of 
all SIX (secreted in xylem) effector genes (Schmidt et al. 2013). 
Similarly members of the M. oryzae Avr-Pita effector gene family are 
flanked by a retrotransposon, presumably contributing to its 
translocation across the genome (Chuma et al. 2011). Presence/ 
absence of polymorphism and translocation of effector genes that 
are associated with unstable genomic regions has also been ob- 
served in asexual lineages of M. oryzae (Yoshida et al. 2009; Chuma 
et al. 2011), and it was hypothesized that parasexual recom- 
bination facilitates the exchange of effector genes between these 
lineages (Noguchi et al. 2006; Chuma et al. 2011). Parasexual re- 
combination involves fusion of haploid fungal hyphae, followed 
by the formation of diploid nuclei through hybridization that can 
later return to a haploid state. Mitotic recombination may result in 
nonsexual exchange of genetic material without meiosis. In par- 
ticular cases, the hybridized nuclei are stable, resulting in a hybrid 
species with a doubled genome such as has been suggested for 
V. longisporum (Inderbitzin et al. 2011). Parasexuality has been 
identified in many, mostly asexually propagating fungi during in 
vitro culturing and was shown to contribute to adaptive capabil- 
ities of the model fungus Aspergillus nidulans in vitro (Schoustra 
et al. 2007). Parasexual recombination may also be involved in the 
exchange of effector genes between asexual lineages of V. dahliae. 
Laboratory studies have shown that transgenic V. dahliae strains 
that carry auxotrophic markers are capable of hyphal fusion and 
exchange of genetic material, depending on their vegetative 
compatibility (Hastie 1964, 1973). Likewise, transfer of pathoge- 
nicity chromosomes between asexual lineages of F. oxysporum may 
be mediated by parasexual cycles (Corell 1991; Ma et al. 2010). 
However, the significance of parasexuality in nature, and the ge- 
netic factors that control parasexual compatibility between line- 
ages, are poorly understood, although evidence for parasexual re- 
combination in field populations of M. oryzae has been reported 
(Zeigler et al. 1997; Glass et al. 2000). 

In our study, chromosomal rearrangements were found to be 
frequently associated with retrotransposons, a class of repetitive 
sequence elements that is known to promote genome instability 
through recurrent excision and insertion as well as by mediating 
homologous recombination during cell division (Zolan 1995; 
Lemoine et al. 2005; Maxwell et al. 2011). In organisms without 
a cycle of sexual reproduction, such as V. dahliae, homologous re- 
combination between highly similar copies of transposable ele- 
ments might help to generate genetic diversity by facilitating mi- 
totic crossover, thereby creating new genetic combinations that are 
potentially beneficial. In addition, altered cosuppression, silencing 
of dispersed homologous genomic regions through directed DNA 
methylation, might act on proliferation of transposable elements to 
affect genome stability (Martienssen and Colot 2001; Selker et al. 
2003; Muszewska et al. 2011). Although it is generally believed that 
asexual reproduction limits genetic variation, and consequently 



1278 Genome Research 

www.genome.org 



Verticillium chromosomal reshuffling 



limits adaptive capability of the organism, we provide evidence for 
chromosomal plasticity as facilitator of asexual genome evolution. 

Methods 

SNP identification 

Illumina reads (de Jonge et al. 2012) were mapped onto the Verti- 
cillium dahliae VdLsl7 reference genome using GSNAP (version 
2012-04-16; Wu and Nacu 2010) with default settings and pro- 
cessed by Picard (version 1.57; http://picard.sourceforge.net). SNPs 
were identified by SAMtools mpileup (version; 0.1.18; http:// 
samtools.sourceforge.net) and filtered using vcfutils.pl varFilter 
(-Q20, -DlOO) and a minimum allele frequency of 0.8. SNP sta- 
tistics were determined by VCFtools (version 0.1.10; Danecek 
et al. 201 1) and variant annotation by SnpEff (version 2.1; http:// 
snpeff.sourceforge.net). To specifically identify genes that are sub- 
ject to positive selection, we applied the Nei Gojobori method, where 
rates of iCa and are compared using a Z-test (Nei and Kumar 2000; 
Stukenbrock and Dutheil 2012). Z-values >2326 are considered 
significant with P < 0.01. 

Coverage breadth was calculated for all genes as the percent- 
age of nucleotides with minimum one read aligned using BEDtools 
(version 2.16.1; Quinlan and Hall 2010). Genes were considered 
absent when breadth was <0.1. 

Maximum likelihood phylogenetic analysis was performed in 
MEGA5 (Tamura et al. 2011). For each strain, as input we extracted 
for all nonredundant SNP positions the respective base call, 
resulting in 11 sequences, each containing 236,785 bp. 

Assembly, alignment, and identification of rearrangements 

Draft assemblies generated previously (de Jonge et al. 2012) were 
used for all strains, except for JR2. Mate-pair library preparation 
and sequencing (50 bp, ~5-kb insert) of JR2 was performed by 
the BGI (Hong Kong). For de novo assembly, Velvet (version 
1.2.03; Zerbino and Birney 2008) was used with cov_cutoff=6, 
exp_cov=auto, k-mer=31 and shortMatePaired=yes to generate 
a new assembly based on the 500-bp paired-end and 5-kb mate-pair 
libraries. Optical mapping, i.e., the construction of ordered ge- 
nome-wide, high-resolution (>150x) restriction maps, was per- 
formed by BGI using the Argus System, and contigs were placed 
using MapSolver (version 3.2; OpGen). 

Whole-genome alignments and dot plots were generated by 
MUMmer 3.0 (Kurtz et al. 2004) using NUCmer with default set- 
tings (except for -115 and-maxmatch) and Mummerplot. We used 
custom scripts to identify rearrangements and associated break- 
points. Core and LS genomic regions were determined by assessing 
coverage breadth of the alignments on 1-kb nonoverlapping 
windows using BEDTools. In addition, we determined the per- 
centage of gaps for each region. Regions were considered lineage 
specific with an alignment breadth of <0.2. Core genomic regions 
were defined based on their minimum presence in all strains ex- 
cept one in order to compensate for assembly and alignment in- 
consistencies. Similarly, genes were considered as core genes when 
identified minimally in all strains except one with an alignment 
breadth of >0.8 for each strain. Remaining genes are considered 
lineage specific. Collinear blocks, repetitiveness, percentage GC, 
and LS genomic regions in the VdLsl7 and JR2 genomes were vi- 
sualized using Circos (version 0.55; Krzywinski et al. 2009) and 
GBrowse (version 2.49; Stein et al. 2002). 

Multiple chromosomal rearrangements were verified by PGR. 
To this end, primer pairs spanning predicted breakpoints were 
designed for JR2 and VdLsl7 to selectively amplify breakpoint re- 
gions in either of the two genotypes and used in PGR reactions on 



genomic DNA. Quality of input DNA and integrity of syntenic 
sequences flanking the breakpoints in JR2 and VdLsl7 were veri- 
fied by control primer sets (Supplemental Table 5). 

Repeat identification 

Repetitive elements were identified and classified by RepeatMod- 
eler (version 1.0.5; http://www.repeatmasker.org/RepeatModeler. 
html), incorporating the two de novo repeat-finding programs 
RECON (version 1.0.7; Bao and Eddy 2002) and RepeatScout 
(version 1.0.3; Price et al. 2005), the tandem repeat finder TRF 
(version 3.2.1), and the annotated repeat library from RepBase (re- 
lease 16.12). Repetitive sequence families were used by RepeatMasker 
(version open-3.3.0; http://www.repeatmasker.org) to identify and 
mask all repetitive sequences applying the sensitive mode. Full- 
length long terminal repeat (LTR) retrotransposons were identified 
by LTR_F1NDER (version 1.0.5; Xu and Wang 2007). Statistical 
significant spatial correlations between repetitive sequences, break- 
points, and genomic regions were assessed by GenometriCorr 
(version 1.1.7; http://genometricorr.sourceforge.net). 

Karyotyping 

V. dahliae mycelium was prepared for protoplasting following the 
mycelium-based fungal biomass preparation method (Mehrabi 
et al. 2012). Mycelium was digested in 1 M sorbitol containing 1% 
(w/v) glucanex and 0.5% driselase (Sigma-Aldrich) at 32°C until 
protoplast concentration reached 10^/mL. Protoplast plugs were 
prepared and stored as described (Mehrabi et al. 2012). Karyotyp- 
ing was carried out using a CHEF Mapper XA pulsed field electro- 
phoresis system (Bio-Rad) using the auto algorithm function with 
low- and high-molecular weight settings at 2 and 6 Mb, respec- 
tively. Chromosomes were separated in 0.8% low EEO agarose (US 
Biological) gels. Chromosome size markers from Hansenula wingei 
and Schizosaccharomyces pombe (Bio-Rad) were included as refer- 
ence. After electrophoresis, gels were stained with ethidium bro- 
mide (1 fxg/mL) in water for 1 h and destained in water for 2 h. 

Gene prediction, expression, and annotation 

Gene predictions on the genome sequences of all V. dahliae strains 
were performed using a combinatorial approach, applying the 
EVidenceModeler (version r2012-06-25; Haas et al. 2008) to com- 
bine protein-coding gene evidence from V. dahliae-tTdined Augustus 
(version 2.6; Stanke et al. 2006, 2008), GeneMarkES (version 2.3c; 
Borodovsky and Lomsadze 2011), V. dahliae-tiained GlimmerHMM 
(version 3.0.1; Majoros et al. 2004), protein alignments to the latest 
Fusarium graminearum (release 3.2, MIPS, Munich) (Wong et al. 
2011) and Magnaporthe oryzae (release 8, Broad Institute) proteins 
using the Analysis and Annotation Tool (AAT) package (version 
r0305201 1; Huang et al. 1997). For deep transcriptome sequencing 
and mapping, ~2 Gb of reads from V. dahliae strain JR2-infected 
Nicotiana benthamiana plants and ~ 1 Gb of reads from V. dahliae 
JR2 cultured in vitro on Czapek Dox medium were mapped on the 
V. dahliae ]R2 genome using TopHat (version 2.0.8; Trapnell et al. 
2009; de Jonge et al. 2012). Cufflinks (version 2.0.2; Trapnell and 
Salzberg 2010) was used to assemble novel transcripts and isoforms 
from the mapped reads that were subsequently incorporated in the 
final gene-set for JR2 using the EVidenceModeler after manual 
curation. Relative expression for each gene in each experiment was 
determined by assessing the number of reads mapping to each 
gene using HTSeq (version 0.5.4pl; www-huber.embl.de/users/ 
anders/HTSeq/) applying the default "union" setting and subse- 
quently reported in the number of Reads Per Kilobase of transcript 
per Million mapped reads (RPKM). For differential gene expression 
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analysis we compared in planta and in vitro expression using Excel 
and CummeRbund (http://compbio.mit.edu/cummeRbund/). The 
predicted proteins were mined for candidate effectors (de Jonge 
2012) by functional annotation using BLASTp analyses against the 
nonredundant database at NCBI (http://www.ncbi.nih.nlm.gov) 
and Pfam domain scanning (Finn et al. 2008) using the Blast2GO 
suite (Conesa et al. 2005). Secreted proteins were predicted by 
SignalP4 (Petersen et al. 2011). Gene density was determined 
as previously described (Raffaele et al. 2010). In short, 5'- and 
3 '-flanking intergenic regions were calculated, scored in two-di- 
mensional bins, and plotted. Filled contour plots were generated 
using R (http://www.r-project.org/). 

Functional analysis 

Knockouts were generated by amplifying the sequences flanking 
coding sequences using the primer sets KO-XLOC_009059-F1/KO- 
XLOC_009059-R1 and KO-XLOC_009059-F2/KO-XLOC_009059- 
R2 fox XLOC_009059, KO-XLOC_008951-F1/KO-XLOC_008951-R1 
and KO-XLOC_00895 1-F2/KO-XLOC_00895 1-R2 for XLOC_008951 , 
and KO-VDAG_05180-F1/KO-VDAG_05180-R1 and KO-VDAG_ 
05180-F2/KO-VDAG_05180-R2 for VDAG_05180 (Supplemental 
Table 4). Likewise, knockouts of XLOC_0090 5 9 dnd XLOC_008951 
were generated by amplification of the flanking sequences (Sup- 
plemental Table 4). PGR products were subsequently cloned into 
pRF-HU2 (Frandsen et al. 2008). V. dahliae transformation and 
subsequent inoculations on tomato (cv. Motelle and MoneyMaker) 
plants to assess the impact on virulence were performed as de- 
scribed (Fradin et al. 2009). Plants were regularly inspected during 
a 2-wk interval and photographed at 8- and 12-d post -inoculation 
(DPI). For biomass quantification, the roots and stem below coty- 
ledons of three plants per V. dahliae genotype were flash-frozen in 
liquid nitrogen. To determine canopy surface areas, we measured 
the leaf surface areas of multiple plants from photographs taken 
from the top using ImageJ. The samples were ground to powder, of 
which an aliquot was used for DNA isolation (Fulton et al. 1995). 
Real-time PGR was conducted with primer sets SlRub-Fl/SlRub-F2 
for tomato RuBisCo and VdGAPDH-F/VdGAPDH-R for V. dahliae 
GAPDH (Supplemental Table 4). For expression analyses, 3-wk-old 
Nicotiana henthamiana plants were inoculated with strain VdLsl7 
as previously described (Fradin et al. 2009), harvested at 4, 8, 12, 
and 16 DPI, and flash-frozen in liquid nitrogen. Total RNA was 
extracted using the RNeasy Kit (Qiagen) and cDNA was synthesized 
by Superscript III (Invitrogen). Real-time PGR was conducted with 
primer sets VdGAPDH-F/VdGAPDH-R for V: dahliae GAPDH and 
qVDAG_05180-Fl/qVDAG_05180-Rl for V. dahliae VDAG_05180 
(Supplemental Table 4). 

Data access 

All genome and transcriptome raw read data have been deposited 
at the NGBI Sequence Read Archive (SRA) (http://www.ncbi.nlm. 
nih.gov/sra) under BioProject PRJNA169154 (http://www.ncbi.nlm. 
nih.gov/bioproject). The V. dahliae strain JR2 draft genome sequence 
is deposited at the NGBI WGA archive database (http://www.ncbi. 
nlm.nih.gov/genome) under accession number APFDOOOOOOOO and 
described under BioProject PRJNA175765. 
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